program VMD !  works with df   based on code drivenDf Runge K.f95 which is in frame of source (phase velocity=0)   

  use fft_wrap, only   : fft_init
  use efield
  use evolve

  use dimensions, only : jx, jz, jv, jDirection

  use parameters, only : restart, nt, dt, iDiag, iData, vPhase, runname

  implicit none

  !-- major arrays ---------------------


  real       :: dfu(jv, jz, jy, jx, jDirection) ! distribution function
  real       :: uxMat(jz, jy, jx, jDirection)   ! velocity, x-component
  real       :: uyMat(jz, jy, jx, jDirection)   ! velocity, y-component
  real       :: uzMat(jz, jy, jx, jDirection)   ! velocity, z-component
  real       :: Ex2D(jz, jy, jx)                ! electric field, x-component 
  real       :: Ey2D(jz, jy, jx)                ! electric field, y-component 
  real       :: Ez2D(jz, jy, jx)                ! electric field, z-component


  !-- suplemental arrays --------------

  real       :: f0vec(jv)                       ! background distribution, gaussian
  complex    :: fvc0T(jv)                       ! FT of f0vec

  real       :: ExFilter(jx)
  real       :: EyFilter(jy)

  !-- variables  ----------------------- 

  real       :: t                               ! time 
  real       :: Courant_in_x                    ! advection Courant number 
  real       :: Courant_in_y                    ! advection Courant number 
  real       :: Courant_in_z                    ! advection Courant number 
  real       :: OmegaBounce

  integer    :: it                              ! iteration count
  integer    :: iargc                           ! argument count


 
!-------------------------------------------------------------------------------------

      !--  reading the run_name from the command line

      if(iargc().eq.0) then
         call getarg(0,runname)
         write(*,*)      'Usage: ',trim(runname),' <run name>'
         stop
      end if
      call getarg(1,runname)


      !-- set parameters, coordinate space, and the filter --

      call fft_init()

      call read_input()

      call set_space_vectors(ExFilter, EyFilter)                       ! output

      call write_derived()

      call evolve_init()
      
      call efield_init()



      !-- restore from restart file or create initial conditions --

      call set_ic( dfu, uxMat, uyMat, f0vec, fvc0T, it)                ! output 


      call efield_compute_E(dfu, ExFilter, EyFilter, &                 ! input 
                            Ez2D, Ex2D, Ey2D, OmegaBounce)             ! output

      t  = it*dt


      !-- write initial diagnostics and binaries --

      if (restart .eq. 0) then

         Courant_in_z = vPhase*dt/dz

         call diagnostics( dfu, Ez2D, Ex2D, Ey2D, f0vec, &             ! inputs 
                           t, Courant_in_z, OmegaBounce)               ! inputs

         call write_dispersion()
         call write_f0_initial(f0vec, dfu)                             ! input 

         call write_bin(dfu, uxMat, uyMat, it, it/iData)               ! input
         call write_bin_extra(dfu, Ex2D, Ey2D, Ez2D, it, it/iData)     ! input

      end if



!-------------------------------------------------------------------------------------
! BEGIN EVOLUTION
!-------------------------------------------------------------------------------------

 
do while (it .le. nt) 

   call evolve_F_advect( dfu, &                               ! input/output
                         uxMat, uyMat, f0vec)                 ! input

   call efield_compute_E( dfu, ExFilter, EyFilter, &          ! input 
                         Ez2D, Ex2D, Ey2D, OmegaBounce)       ! output

   call evolve_F_accel(  dfu, &                               ! input/output
                         Ez2D, fvc0T)                         ! input

   call evolve_Uz(       dfu, f0vec, &                        ! input
                         uzMat)                               ! output

   call evolve_U_advect( uxMat, uyMat, &                      ! input/output
                         uzMat)                               ! input
 
   call evolve_U_accel(  uxMat, uyMat, &                      ! input/output
                         Ex2D, Ey2D)                          ! input

   call evolve_U_filter(  uxMat, uyMat, &                     ! input/output
                          ExFilter, EyFilter)                 ! input

   call efield_evolve_E0()                                    ! input


   t  = t + dt
   it = it + 1
   
   Courant_in_z = maxval(abs(uzMat))*dt/dz
   Courant_in_x = maxval(abs(uxMat))*dt/dx
   Courant_in_y = maxval(abs(uyMat))*dt/dy

   !-- diagnostics --

   if ( (it/iDiag)*iDiag.eq.it ) then      

      call diagnostics( dfu, Ez2D, Ex2D, Ey2D, f0vec, &        ! inputs 
                        t, Courant_in_z, OmegaBounce)          ! inputs

   end if

   !-- output binary data for restart and post-processing --

   if ( (it/iData)*iData.eq.it ) then      
     call write_bin(dfu, uxMat, uyMat, it, it/iData)           ! input
     call write_bin_extra(dfu, Ex2D, Ey2D, Ez2D, it, it/iData) ! input
   end if

   !-- check Courant conditions --
      
   if(  Courant_in_z > 1.0 ) then
      write(*,*)   'Courant for z advection =', Courant_in_z, '  is > 1, stop evolution'
      go to 314
   endif
            
   if(  Courant_in_x > 1.0 ) then
      write(*,*)   'Courant for x advection =', Courant_in_x, '  is > 1, stop evolution'
      go to 314
   endif

   if(  Courant_in_y > 1.0 ) then
      write(*,*)   'Courant for x advection =', Courant_in_y, '  is > 1, stop evolution'
      go to 314
   endif

end do      ! END TEMPORAL LOOP
 
314 continue


    write(*,*) 'Dvvvv should be > E_max / p_max^3 = ', maxval(Ez2D) / pvec(np/2)**3
    write(*,*) 'end temporal evolution'  


    !-- write some specrtral spectrum quantities --

    call efield_write_spectra()

    call write_umat_spectra(uzMat, uxMat, uyMat)      ! input

    call write_density_spectra(dfu, fvc0T)            ! input
       
    stop;

end


!-------------------------------------------------------------------------------------
! subroutine deAlias
!-------------------------------------------------------------------------------------

      
subroutine deAlias(jz,et); implicit none;  integer i,jz;  complex et(jz)
if( jz<6 ) return
       do i=jz/2  - jz/6, jz/2 + 2 + jz/6;   et(i)=(0.0,0.0);  end do
return;  end



!-------------------------------------------------------------------------------------
subroutine write_bin(dfu, uxMat, uyMat, it, filenum)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only  :  runname

  implicit none

  real,    intent(in)    :: dfu(jv, jz, jy, jx, jDirection)
  real,    intent(in)    :: uxMat(jz, jy, jx, jDirection)
  real,    intent(in)    :: uyMat(jz, jy, jx, jDirection)
  integer, intent(in)    :: it
  integer, intent(in)    :: filenum


  character*4 :: file_ext

  file_ext = '0000'
  write(file_ext,"(i4.4)") filenum

  open(91, file= trim(runname)//'.dfu.'//file_ext, &
         form='unformatted', access='stream')
  write(91) dfu
  close(91)

  open(92, file= trim(runname)//'.ux.'//file_ext, &
       form='unformatted', access='stream')
  write(92) uxMat
  close(92)

  open(92, file= trim(runname)//'.uy.'//file_ext, &
       form='unformatted', access='stream')
  write(92) uyMat
  close(92)

!  open(93, file= trim(runname)//'.tag')
!  write(93,*) it, filenum 
!  close(93)
 

 end subroutine write_bin


!-------------------------------------------------------------------------------------
subroutine read_bin(dfu, uxMat, uyMat, it)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only  :  runname

  implicit none

  real, intent(out)       :: dfu(jv, jz, jy, jx, jDirection)
  real, intent(out)       :: uxMat(jz, jy, jx, jDirection)
  real, intent(out)       :: uyMat(jz, jy, jx, jDirection)

  integer, intent(out)    :: it
  integer                 :: filenum

  character*4 :: file_ext

  open(93, file= trim(runname)//'.tag')
  read(93, *) it, filenum 
  close(93)

  write(*,*) it, filenum


  file_ext = '0000'
  write(file_ext,"(i4.4)") filenum

  open(91, file= trim(runname)//'.dfu.'//file_ext, &
       form='unformatted', access='stream')
  read(91) dfu
  close(91)

  open(92, file= trim(runname)//'.ux.'//file_ext, &
       form='unformatted', access='stream')
  read(92) uxMat
  close(92)

  open(92, file= trim(runname)//'.uy.'//file_ext, &
       form='unformatted', access='stream')
  read(92) uyMat
  close(92)


end subroutine read_bin

!-------------------------------------------------------------------------------------
subroutine write_bin_extra(dfu, Ex2D, Ey2D, Ez2D, it, filenum)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only  :  runname

  implicit none

  real,    intent(in)    :: dfu(jv, jz, jy, jx, jDirection)
  real,    intent(in)    :: Ex2D(jz, jy, jx)
  real,    intent(in)    :: Ey2D(jz, jy, jx)
  real,    intent(in)    :: Ez2D(jz, jy, jx)

  integer, intent(in)    :: it
  integer, intent(in)    :: filenum

  real                   :: dens(jz, jy, jx)

  character*4            :: file_ext

  integer                :: ix, iy, iz


  do ix=1,jx
     do iy=1,jy
        do iz=1,jz
           dens(iz,iy,ix) = sum(dfu(:,iz,iy,ix,:))*dv 
        end do
     end do
  end do



  file_ext = '0000'
  write(file_ext,"(i4.4)") filenum

  open(91, file= trim(runname)//'.Ex.'//file_ext, &
         form='unformatted', access='stream')
  write(91) Ex2D
  close(91)

  open(92, file= trim(runname)//'.Ey.'//file_ext, &
       form='unformatted', access='stream')
  write(92) Ey2D
  close(92)

  open(92, file= trim(runname)//'.Ez.'//file_ext, &
       form='unformatted', access='stream')
  write(92) Ez2D
  close(92)
 
  open(92, file= trim(runname)//'.dn.'//file_ext, &
       form='unformatted', access='stream')
  write(92) dens
  close(92)

 end subroutine write_bin_extra
!-------------------------------------------------------------------------------------
